1. IntroductionOne of the biggest challenges and unsolved problems in condensed matter and quantum field theory is the study of quantum many-body systems, whose state space is exponentially large. This has severely delayed our understanding of many fascinating strongly correlated phenomena, including high temperature superconductivity and quantum spin liquids. Quantum Monte Carlo (QMC) simulation is one of the most successful methods in simulating quantum many-body systems, but fails in the simulation of interacting fermions and frustrated spin models due to the infamous minus sign problem. In recent years, tremendous progress has been achieved in the development of numerical renormalization-group methods based on tensor-network states (TNS),[1–11] which have emerged as a powerful theoretical tool for investigating low-dimensional quantum lattice models.
The TNS formulation is a variational ansatz for the ground state. It reduces the dimension of the Hilbert space, which grows exponentially with the system size, to a polynomial growth. Commonly used TNS include the one-dimensional matrix product state (MPS),[12] which is a class of states underlying the density-matrix renormalization group (DMRG),[13] and the two-dimensional projected entangled pair state (PEPS).[3] The accuracy of TNS is controlled by the virtual bond dimension of the local tensors, D. The larger the bond dimension is, the more accurate a TNS is. However, the cost for computing a TNS, especially a PEPS or a projected entangled simplex state (PESS),[11] rises rapidly with increasing D. For example, the minimal cost scales as D12 for PEPS. This has limited the bond dimension that can be handled to be generally less than 13 in two dimensions.[6,10,11] Furthermore, although both MPS and PEPS satisfy the area law of entanglement entropy,[14] in a critical or interacting fermion system with a finite Fermi surface, there is a logarithmic correction to the entropy. To describe correctly this logarithmic behavior, a more complex TNS structure is required, namely, the multi-scale entanglement renormalization ansatz (MERA) in one dimension[5] or the branching MERA in two dimensions.[15] The cost for handling these MERA-type wavefunctions is even higher. Resolving this difficulty requires a new approach that can improve significantly the accuracy of TNS without relying on the increase of the bond dimension.
In this work, we propose a generalized Lanczos method to solve quantum lattice models using TNS. This method is an adaptation of the Lanczos method for the tensor network algorithm, which generates a set of orthonormal many-body basis states (i.e., the Krylov subspace), represented using TNS, by applying the Hamiltonian to the iteratively generated basis states. At each iteration step, a new TNS is generated by minimizing a cost function. However, as a TNS is only an approximate representation of a quantum state, the Hamiltonian is not tri-diagonalized unlike in the standard Lanczos method. By diagonalizing the Hamiltonian in this set of basis states, a better ground state is obtained and represented as a linear superposition of all the generated TNS.
This paper is organized as follows. We first introduce the TNS-Lanczos method in Section 2. In this section, we also present a detailed discussion on how to minimize a cost function to obtain a Lanczos vector represented by a TNS. In Section 3, we present the results obtained with the TNS-Lanczos method for the antiferromagnetic Heisenberg model in both one and two dimensions. A summary is given in Section 4.
2. MethodThe TNS-Lanczos method starts from a TNS, |ψ1⟩, that is determined by variationally minimizing the energy functional
subject to the constraint ⟨
ψ1|
ψ1⟩ = 1. |
ψ1⟩ is set as the first normalized many-body basis state, |
Ψ1⟩ = |
ψ1⟩. A new TNS, |
ψ2⟩, with the same bond dimension is then generated by minimizing the cost function
where
h11 is the energy of the basis state |
Ψ1⟩, defined by Eq. (
6). From |
ψ2⟩, we can construct a new normalized basis state |
Ψ2⟩ that is orthogonal to |
Ψ1⟩,
where
a2 is the normalization constant.
Similarly, from |Ψ1⟩ and |Ψ2⟩, we can generate another TNS, |ψ3⟩, by minimizing a cost function similar to Eq. (2). Again, from |ψ3⟩, we can construct a new normalized basis state |Ψ3⟩ that is orthogonal to both |Ψ1⟩ and |Ψ2⟩. By continuing this iteration, we create a set of orthonormal basis states {|Ψα⟩; α = 1, . . ., k} with ⟨Ψα|Ψβ⟩ = δα,β.
In general, to find the basis state |Ψα+1⟩, we first generate a TNS, |ψα+1⟩, by minimizing the cost function
where
hαβ is the matrix element of the Hamiltonian in this set of basis states,
and
|
Ψα+1⟩ is defined by
with
aα+1 being the normalization constant.
The cost function defined in Eq. (2) or more generally in Eq. (4) can be minimized iteratively by sweeping over all local tensors, similar to in the full-update calculation.[7] At each step, a local tensor Ai at site i is updated with all the other tensors fixed. The corresponding cost function for this local tensor can be expressed as a quadratic function of Ai
where
Ni is a matrix determined by contracting the inner product ⟨
ψα+1|
ψα+1⟩ excluding
Ai and its Hermitian conjugate
, and
bi is a vector determined by ⟨
ψα+1|Φ
α⟩ excluding
Ai. Minimizing this function with respect to
Ai yields a set of linear equations
Ai is determined by solving this set of equations. Repeating this iterative step until the cost function has converged, we then obtain |
ψα+1⟩.
|Φα⟩ is a linear superposition of α TNS. The cost for contracting ⟨ψα+1|Φα⟩ scales linearly with α. Thus to generate n orthonormal TNS basis states, the total computational time needed increases roughly quadratically with n. But the contraction of ⟨ψα+1|Φα⟩ can be readily parallelized, which can reduce the total computational time from n2 to n.
From the above iteration, we find k TNS {|ψ1⟩, |ψ2⟩, . . ., ψk⟩} and k orthonormal basis states {|Ψ1⟩, Ψ2⟩, . . ., Ψk⟩} in one round of Lanczos iterations. In this basis space, the Hamiltonian can be represented as a k × k matrix whose matrix elements are defined by Eq. (6). By diagonalizing this matrix, we obtain the ground-state energy and the corresponding eigenfunction. However, the computational time scales as (k − 1)2, and furthermore, due to the accumulation of machine errors, the orthogonality of the basis states so generated may be lost when k becomes large. To avoid these problems, it is better not to use a large k in the calculation, and instead, we use a relatively small k but repeat the above steps many times. Each time we set the ground state obtained from the previous cycle as the initial basis state, |Ψ1⟩. Of course, starting from the second round of iterations, the initial basis state is no longer a single TNS. Instead it is a linear superposition of all the TNS obtained previously. This restarted Lanczos iteration can be repeated many times until the ground-state energy converges. If we use n to denote the number of restarted Lanczos iterations, then the total number of TNS generated in this way equals N = n(k − 1) + 1. The ground-state wave function is a linear superposition of these N TNS, {|ψα⟩, α = 1, . . ., N}, which we call a Lanczos-generated TNS (abbreviated as LTNS).
A LTNS can be generally expressed as
where
cα is the coefficient of the ground-state wave function in this tensor-network representation. This wave function can also be represented as a single TNS. If
is the local tensor of
cα |
ψα⟩, with
mi the quantum number of the physical basis state at site
i, then the local tensor of |
Ψ⟩ at the same site
i,
Ai[
mi], is simply a block-diagonal tensor defined by
The bond dimension of
Ai equals
ND. Clearly, |
Ψ⟩ contains more variational parameters than each pure TNS |
ψα⟩. Thus it is not surprising that it can be more accurate than the wave function |
ψ1⟩ obtained by simply minimizing the ground-state energy.
The memory space needed for storing a LTNS scales linearly with N. The computational time required for generating these basis states scales as n2(k − 1)2 = (N − 1)2. The converged ground-state energy in the large-N limit depends on the total number of TNS generated, but does not depend much on the value of k. In the case that only the ground state is studied, it is sufficient to take k = 2. For a larger k, the ground state energy can converge faster than the k = 2 case during the first tens of iterations, but the entire cost is higher. The results presented below are all obtained with k = 2.
3. ResultsWe test the method using the spin-1/2 antiferromagnetic Heisenberg model with open boundary conditions in both one and two dimensions. The Hamiltonian of the Heisenberg model reads
where
Si is the SU(2) spin operator defined at site
i and ⟨
ij⟩ represents summation over nearest-neighbor lattice sites.
3.1. MPS-based Lanczos calculationsWe first carry out the calculation based on the MPS representation of the basis states. Figure 1 shows how the ground state energy E(n) and the entanglement entropy S(n) vary with the restarted Lanczos number n for the Heisenberg model, obtained with MPS. As the initial E(n) and S(n), i.e. E(0) and S(0), are just the results obtained by the standard DMRG method, the improvement of our method over the DMRG is quite significant if the same bond dimensions are used.
The entanglement entropy grows very rapidly with n in the first tens of iterations and converges to a constant in the limit n → ∞. The ground-state energy is strongly anti-correlated with the entanglement entropy, dropping quickly with increasing n. For example, the relative error in the ground-state energy for L = 20 is reduced by nearly two orders of magnitude for D = 8 and three orders of magnitude for D = 12 at n = 300. The ground-state energy keeps descending with increasing n, but with a smaller slope when the entanglement entropy becomes saturated.
By directly comparing the entanglement entropy of the initial MPS |Ψ1⟩ to that of the converged LTNS, shown in Fig. 2(a), we find that the latter contains much more entanglement than the former. In particular, the entanglement entropy of the LTNS with D = 8 already reaches the values calculated with DMRG by keeping as many as 1024 states, which can be regarded as quasi-exact. Furthermore, the entanglement entropy of the LTNS varies logarithmically with the system size, indicating that the LTNS can describe correctly the scaling behavior of a critical system, even though each individual MPS is bounded by the entanglement area law.
Figure 2(a) shows the entanglement entropy as a function of the lattice size for the Heisenberg spin chain. For all the lattice sizes we have studied, we find that the converged entanglement entropy agrees with the equation derived from the conformal field theory[18]
where
b is a non-universal constant and
c is the central charge. In our calculation, both
c and
b are weakly dependent on the bond dimension. For
D = 8 and 12,
c is found to be approximately 1.014, close to the exact result,
c = 1.
The entanglement entropy of the LTNS, as revealed in Fig. 1(a), is contributed primarily by the first few tens of TNS generated by the Lanczos iterations. To reproduce the logarithmic L-dependence, the number of these TNS that have the most significant contribution to the entanglement entropy should increase with the lattice size. This is indeed what we find. Figure 2(b) shows how the entanglement entropy S(n) approaches its converged value S(∞), with n for several different values of L. Clearly, the number of TNS that have the most significant contributions to the entanglement entropy increases with L. If n0 is the number of TNS whose contribution to the entanglement entropy is within 1% of the converged value, i.e., 1 − S(n0)/S(∞) ≤ 1%, we find that n0 varies almost linearly with the system size (inset in Fig. 2(b)).
The correlation function is an important indicator of low energy properties of the ground state. From its long-range behavior, one can judge whether the system is gapped or gapless. However, the long-range correlation function is generally difficult to determine accurately. Figure 3 compares the correlation function and the entanglement entropy obtained using our method with those obtained by the DMRG calculations for the S = 1/2 Heisenberg spin chain. We calculate the correlation function when the entanglement entropy obtained using the TNS-Lanczos with D = 8 becomes saturated at n = 400. Our results agree very well with those obtained from the DMRG calculation by keeping D = 1024 states. This indicates that our method can yield accurate results not just for the ground state energy and the entanglement entropy, but also for the long range correlation functions with relatively very small bond dimensions.
We have also calculated the ground state energy E(n) as a function of n for the two-dimensional Heisenberg model on the 10 × 10 lattice using the MPS-based Lanczos method. The results are shown in Fig. 4. By keeping D = 1000 states, we find that the ground state energy obtained by taking n ∼ 50 iterations is already lower than the DMRG result obtained with D = 2000. For this system, D = 4000 is close to the upper limit of the bond dimension in the standard DMRG calculation without implementing the SU(2) symmetry due to the limitation of memory space. In this case, we can still improve the result by taking a number of Lanczos iterations. Thus, as revealed by the figure, our method can be used to effectively enlarge the maximal bond dimension of MPS with very limited increase of computational cost.
3.2. PEPS-based Lanczos calculationsFor the ground state represented using the PEPS or other two-dimensional TNS, the improvement of our method over the conventional tensor-network algorithm is even more pronounced. Figure 5 shows how the ground-state energy and the entanglement entropy obtained with PEPS vary with the Lanczos number for the isotropic Heisenberg model on the 4 × 4 square lattice. Similar to the one dimension case, the ground-state energy drops very rapidly in the first tens of iterations while the entanglement entropy grows most rapidly. As the ground state of the two-dimensional Heisenberg model satisfies the entanglement area law, we find that the ground-state energy E(n) converges even faster than its one-dimensional counterpart. For example, with just 200 Lanczos iterations, the accuracy of the ground-state energy is improved by 9 orders of magnitude for D = 4. Moreover, the ground state energy converges exponentially to the exact result Eex in the large n limit
where
a and
μ are two parameters depending on both the bond dimension and the system size.
μ measures the rate of convergence. Given
L,
μ increases rapidly with increasing
D. For
D = 4,
μ is quite close to that of the Lanczos ED method.
The TNS-Lanczos method also works very well in a system with larger lattice size or with larger bond dimensions, in which the overlap between two PEPS wave functions is calculated approximately using the boundary MPS or TEBD method.[7,19,20] As shown in Fig. 6, the ground state energies obtained with our method for the Heisenberg model on the 10 × 10 lattice in the PEPS representation with D = 2 and 3 are already lower than those obtained by the full-update method with D = 3 and D = 6, respectively. The largest bond dimension χ of the boundary MPS used in our calculation is 80, which is large enough to obtain converged physical expectation values. In Fig. 7 we show the convergency of the ground state energy with χ. More results obtained with our method on other systems sizes are shown in Fig. 8.
Figure 9 shows a power-law fitting for the ground state energy per site obtained by the PEPS-based Lanczos method with D = 3 and n = 50,
The extrapolated ground state energy in the thermodynamic limit is
E0 = −0.66812(74), which is slightly lower than the value obtained by the SU(2)-invariant PEPS calculation on the infinite lattice
[21] with
D = 7, −0.6677, but slightly higher than the result of Monte Carlo simulation,
[22] −0.6694421(4). The exponent
α is found to be 0.91717, close to but slightly lower than 1, which is the leading finite-size exponent expected for an open boundary system.
4. SummaryOur proposed generalized Lanczos method provides a powerful numerical tool to solve quantum lattice models using TNS. It improves significantly the existing tensor-network algorithm and allows the ground state to be more accurately calculated using TNS. At each step of Lanczos iteration, the bond dimension of each newly added TNS is unchanged, but the number of parameters is increased, which allows us to obtain a lower, hence better, ground state energy. This implies that the TNS-Lanczos method can effectively enlarge the maximal bond dimension of TNS that can be handled, especially for PEPS or other TNS in two dimensions. Moreover, the ground state wave function obtained with this method contains more entanglement than a single TNS. It can describe correctly the logarithmic correction to the area law of entanglement entropy in a critical system without invoking a mulitscale entangled TNS, such as MERA.
In a DMRG calculation (similarly in other TNS-based calculations), there is always an upper bound on the maximal bond dimension or the number of states retained, Dmax, that can be handled due to the fast growth of both the computational time and the memory space with Dmax. However, we can effectively enlarge the maximal bond dimension by taking a number of TNS-Lanczos iterations, starting from the DMRG ground state wave function by keeping Dmax basis states. The cost for carrying out this Lanczos calculation is small, because it can be controlled just to grow linearly with the number of iterations. It indicates that in a DMRG or other TNS calculation where the bond dimension already reaches its maximal value, we can still significantly improve the results by taking the Lanczos iterations. This is an advantage of this method in comparison with other TNS methods, especially in the calculation of two-dimensional quantum lattice models with PEPS or PESS where the bond dimension that can be treated is very small.
In Refs. [23] and [24], a set of Lanczos-generated MPS were used to compute dynamic correlation functions in one dimension. In those works, each quantum many-body basis state is approximately represented by a single MPS that is
determined simply by using the standard Lanczos tridiagonal formulism. More explicitly, |ψα+1⟩ is obtained by truncating a higher-dimensional MPS constructed from (H − hαα)|ψα⟩ and |ψα−1⟩, rather than by minimizing the cost function defined by Eq. (4). The basis states such generated, as noted in Ref. [24], suffer severely from the non-orthogonality problem, which may cause a large error in the final result. This problem can be remedied by using our approach.
The TNS-Lanczos approach can in principle be applied to the MPS, PEPS, or other kind of TNS with any kind of boundary conditions. It can be extended to calculate the second or even higher excitation states and the energy gap by targeting two or more basis states at each Lanczos iteration. This can be regarded as a generalization of the block Lanczos method. It can be extended to finite temperature[25] and to a non-Hamiltonian system to compute, for example, thermodynamic quantities using quantum transfer matrices.[26,27] Other kind of Krylov subspace methods similar to the Lanczos method, for example, the Arnoldi[28] and the conjugate-gradient[29] methods, can also be used to generate the Krylov basis states.
